Repository navigation
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #411 +/- ##
==========================================
- Coverage 92.48% 92.33% -0.15%
==========================================
Files 31 31
Lines 3139 3157 +18
==========================================
+ Hits 2903 2915 +12
- Misses 236 242 +6 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
This is how benchmark results would change (along with a 95% confidence interval in relative change) if 503361b is merged into master:
|
jgabry
left a comment
jgabry
left a comment
There was a problem hiding this comment.
Thanks, this looks pretty good. I had Codex check that the k values from sis()/tis() are identical to psis() (example data and a simulated 4000×2000 matrix with one heavy-tailed column) and that the weights are unchanged. I made a few suggestions, mostly about side effects for users.
And one other thing: I think the cost is more than "a little". I asked Codex to check and with S = 4000 and N = 2000, sis() went from 0.38s to 1.66s and tis() from 0.63s to 1.60s, about the same as psis() (1.69s). That's probably fine for someone calling it directly, but loo_subsample(loo_approximation = "tis"/"sis") calls loo.function() with these methods and throws away the k values. With N = 5000, S = 4000 and 200 subsampled observations it went from 5.5s to 7.6s. Do you think that's fine? I'm not sure how much the slowdown would be for bigger examples, maybe a lot? Is it worth skipping the diagnostic in that case?
| k_threshold <- ps_khat_threshold(S) | ||
| tail_len <- n_pareto(r_eff, S) | ||
| # all methods report the Pareto k diagnostic, which needs a long enough tail | ||
| throw_tail_length_warnings(tail_len) |
There was a problem hiding this comment.
With this change, a user calling sis() or tis() users the multi-column tail-length warning, which ends with "Skipping the following columns: 1, 2, ...". But nothing is skipped for these methods except the diagnostic, because the weights are computed as usual, so this could maybe be a confusing warning message? Maybe not a big deal. We could maybe reword the message in throw_tail_length_warnings() to something like e.g. "Pareto k is set to Inf for the following columns: ..."? I'm not sure it's worth changing that message though. Just pointing this out but not a big deal.
| #' * `pareto_k`: Not used in `sis`, all set to 0. | ||
| #' * `pareto_k`: Estimates of the shape parameter \eqn{k} of the | ||
| #' generalized Pareto distribution, computed from the log ratios as in | ||
| #' [psis()]. `sis()` uses them only as a diagnostic. |
There was a problem hiding this comment.
sis() and tis() now warn using the PSIS threshold min(1 - 1/log10(S), 0.7) and pareto_k_table() calls everything below the threshold "good". But for plain IS don't the ratios already have infinite variance when k > 0.5? So e.g. with S = 1000 a k of 0.6 is "good" for sis() and gives no warning. Is the PSIS threshold what we want for sis/tis? Might be simplest to keep it but put in the docs that the threshold and the labels in pareto_k_table() are the PSIS ones? What do you think?
| #' * `pareto_k`: Not used in `tis`, all set to 0. | ||
| #' * `pareto_k`: Estimates of the shape parameter \eqn{k} of the | ||
| #' generalized Pareto distribution, computed from the log ratios before | ||
| #' truncation as in [psis()]. `tis()` uses them only as a diagnostic. |
There was a problem hiding this comment.
Same comment that I made in sis.R.
| cat( | ||
| "Computed from", dims[1], "posterior draws and", | ||
| dims[2], "log-likelihood terms using", class(x)[1], ".\n" | ||
| dims[2], "log-likelihood terms using", paste0(class(x)[1], ".\n") |
There was a problem hiding this comment.
Not about this line, but related: print.importance_sampling_loo() (line 57) still doesn't print the Pareto k table. So loo(x, is_method = "sis") will now warns "Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic')", but the printed output has no diagnostic. If k means something for sis/tis now, should we print pareto_k_table() like print.psis_loo() does? I'm not sure actually, and it doesn't necessarily need to be decided in this PR
| #' @return Scalar Pareto k estimate. `Inf` if the tail is too short or constant, | ||
| #' or if the fit fails, as in `do_psis_i()`. | ||
| #' | ||
| pareto_k_i <- function(log_ratios_i, tail_len_i) { |
There was a problem hiding this comment.
This repeats the tail selection, cutoff and constant-tail check from do_psis_i(). do_psis_i(log_ratios_i, tail_len_i)$pareto_k gives identical k (I had Codex check on 2000 columns) and would remove the duplicated code. Although the diagnostic step would take longer because of the smoothing. I'm not entirely sure which to prioritize, so I'm ok with keeping this as is or remove the duplicate code.
| lw_sorted <- sort.int(log_ratios_i - max(log_ratios_i)) | ||
| tail_ids <- seq(S - tail_len_i + 1, S) | ||
| lw_tail <- lw_sorted[tail_ids] | ||
| if (abs(max(lw_tail) - min(lw_tail)) < .Machine$double.eps / 100) { |
There was a problem hiding this comment.
The constant-tail case isn't covered by the tests. test_psis.R:146 tests it for do_psis_i(), so we could we add the same case for sis()/tis()
|
Based on your comments, I think we should actually drop this and put NA for Prateo k for these |
3 participants
Fixes #227.
This PR was assisted by Claude
Summary
sis()andtis()setpareto_kto 0 for every observation, sopareto_k_values()reported the best possible diagnostic value, although loo had not computed any diagnostic.psis()reportsInfwhen the tail has fewer than 5 draws, so for few draws the three methods disagreed (Infforpsis(), 0 for the other two). With many draws they disagreed too:psis()returned the estimated k, and the other two still returned 0.The Pareto k diagnostic depends only on the log ratios, because
psis()estimates k before it smooths the tail.sis()andtis()now compute the same k aspsis()from the same log ratios. The weights of all three methods do not change.Change
pareto_k_i()inR/psis.Rfits the generalized Pareto distribution to the same tail asdo_psis_i(), with the same cutoff and the sameposterior::gpdfit()call, and skips the smoothing. It returnsInfin the cases wheredo_psis_i()does: a tail shorter than 5 draws, a constant tail, or a failed fit.do_sis_i()anddo_tis_i()usepareto_k_i().tis()computes k from the log ratios before truncation.do_importance_sampling()gives the "Not enough tail samples" warning for all three methods. Before, onlypsis()gave it.sis()andtis()now give the high Pareto k warning when k is above the threshold, aspsis()does.pareto_kandtail_lenforsis()andtis()describes the diagnostic.print()of asis_looortis_looobject wrote "using sis_loo ." with a space before the period. I removed the space.For the example in #227 (3 iterations of 2 chains), all three methods now return
Inffor all 32 observations and give the same warnings. With 500 and 1000 draws ofexample_loglik_array(), the three methods return identical k values.sis()andtis()now fit the generalized Pareto distribution once per observation, aspsis()does, so they take a little more time.Tests
tests/testthat/test_tisis.R:psis, sis and tis report the same Pareto k diagnostic: identical k forexample_loglik_array(). An observation with t-distributed log ratios (1 degree of freedom) gives k > 0.7 and the high-k warning for each method.psis, sis and tis report Inf for too few draws: the Inconsistent Pareto k-values for SIS? #227 example givesInffor all observations and both warnings for each method.sis()andtis()with 16 and 10 draws. They expected no warning; they now expect the tail-length warning, the same onepsis()gives.With the master R files, the new and changed tests give 16 failures. With this PR,
devtools::test()withNOT_CRAN=truepasses all 1283 tests and gives no warnings (R 4.6.1, testthat 3.3.2, Linux x86_64).